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INTRODUCTION 


Recent  work  on  surface  waves  has  shown  that  it  is  possible  to 
detect  weak  surface  wave  signals  at  teleseismic  distances  by  means 
of  a  matched  filter  operation  whereby  the  long  dispersed  waveforms 
are  compressed  into  autocorrelatj.on-liko  waveforms.  In  that  opera¬ 
tion  the  matched  filtered  output  is  the  crosscorrelation  of  a 
reference  signal  with  the  smaller  event  of  interest. 

The  objective  of  the  present  study  is  to  investigate  methods 
for  reliably  estimating  the  spectra  of  weak  surface  waves  in  the 
presence  of  noise.  In  particular  we  attempt  to  take  advantage  or 
the  signal  enhancement  that  can  be  achieved  by  Soth  matched  filter¬ 
ing  and  array  summing.  Reliable  spectral  estimates  are  important 
because  the  spectral  differences  in  surface  wave  excitation  observed 
for  large  explosions  (m^  >  5)  and  earthquakes  of  comparable  body 
wave  magnitude  are  expected  to  exist  for  very  small  events  also. 
Moreover,  reliable  spectral  estimates  are  essential  for  using  the 
frequency  dependence  of  surface  wave  radiation  patterns  as  a  diag¬ 
nostic  (earthquakes  are  expected  to  exhibit  frequency  dependent 
radiation  patterns  while  explosions  should  not). 

In  this  study  we  have  developed  new  methods  of  obtaining 
spectra  using  both  single  channel  and  multichannel  processing.  Tor 
both  test  cases  and  actual  events  we  compare  the  results  using  these 
new  methods  with  conventional  spectral  estimates.  As  expected,  the 
results  are  essentially  the  same  f^or  all  methods  until  low  signal 
to  noise  ratios  are  encounterfd  (S/N  *  2  for  the  test  cases),  where 
the  new  methods  perform  more  reliably. 

Computer  programs  necessary  to  obtain  spectra  for  both  individ¬ 
ual  stations  and  arrays  using  these  methods  have  been  written.  Con¬ 
ventional  and  matched  filter  spectra  are  obtained  individually  and 
for  the  array  sum.  The  program  has  many  options  which  3llow  one  to 


use  an  observed  signal  as  the  reference  filter,  to  generate  a 
theoretical  reference  signal,  to  adjust  a  given  reference  signal 
for  differences  in  epicentral  distance,  or  to  use  a  chirp  filter 
as  the  matched  filter.  In  addition,  event  spectra  may  be  saved  on 
tape  so  that  eventually  region-by-region  comparisons  of  surface 
wave  excitation  as  a  function  of  magnitude,  depth,  etc.  can  be 
obtained  routinely. 

In  the  sections  which  follow  we  develop  the  methods,  analyze 
several  test  cases  and  present  results  for  actual  teleseismic  events 
observed  at  LASA  to  illustrate  the  techniques.  An  analysis  of  surface 
wave  spectra  for  many  events  for  different  regions  utilizing  these 
techniques  will  be  published  in  a  separate  report. 


METHODS 


The  standard  method  of  computation  of  a  power  spectral 
estimate  for  a  signal  immersed  in  noise  is  to  compute  the  smoothed 
power  spectrum  of  the  time  window  containing  the  signal.  From  this 
smoothed  power  spectrum,  a  noise  correction  can  be  subtracted  if 
the  noise  can  be  assumed  to  be  random  and  stationary.  This  noise 
correction  is  obtained  from  the  smoothed  power  spectrum  of  a  noise 
sample  n2(t)  preceding  the  signal.  That  is,  if 

x(t)  =  s  (t)  ♦  n1  (t) 


then 


Pc  (w)  S  Pv(w)  -  (Tv/T  )  P  (w) 

n2  n2 


or 


U) 


Ps  C^)  *  I  X  (•*>)  |  2  -  (Tx/Tn2)  I  N2  (CJ)  I  2 


S2fw)  ♦  2  Re  (S(u>)N*  («))  ♦  iNjCuOl2 


(VTn2)lN2^ 


=  S‘(u)  +  n(w) 
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where 


x(t) 

S2(w) 

V“) 

X(u) 

Px<“> 

Tx 

N2(w) 

V“ 

?n2 

n(w) 


■  time  series  consisting  '  signal,  s(t)  plus 
noise,  (t) 

=  true  power  spectrum  of  s;t) 

=  estimated  power  spectrum  of  s(t) 

*  Fourier  transform  of  x(t) 

=  power  spectrum  of  x(t) 

■  length  of  window  for  x(t) 

=  Fourier  transform  of  signal-free  noise  sample,  n?(t) 

*  power  spectrum  of  n2(t) 

=  length  of  n2  ( t )  window 

=  residual  noise  power  representing  errors  in  the  noise 
correction. 


If,  however,  the  signal  is  dispersed  and  we  have  some  knowledge 
of  the  dispersion,  as  is  the  case  with  teleseismic  Rayleigh  and  Love 
waves,  it  is  possible  to  make  an  alternate  (improved)  estimate  of 
the  signal  spectrum.  By  matched  filtering  the  seismogram  x(t)  with  a 
large  reference  surface  wave  signal  yQ(t)  from  the  same  source  region, 
the  signal  dispersion  is  removed,  i.e.,  the  signal  energy  is  compressed 
into  a  short  autocorrelation- like  pulse  (Alexander  and  Rabenstine,  1967; 
Appendix  A).  The  Fourier  transform  of  a  short  window  of  this 
filter  output  containing  the  signal,  divided^  by  the  amplitude  response 
of  the  reference  filter,  will  give  a  smoothed  spectral  estimate  of 
the  desired  signal  s(t).  Since  the  length  of  the  window  in  the 
matched  filtered  output  required  to  include  the  majority  of  the  signal 
is  significantly  shorter  than  the  original  dispersed  signal  length, 
and  since  the  mean  and  RMS  level  of  stationary  random  noise  is 
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unchanged  by  the  matched  filtering,  the  signal  power  spectrum 
estimated  from  this  window  T  will  have  proportionally  (T  /T  )  less 
noise  contamination  than  that  computed  directly  from  the  signal 
window  Tx  (see  Appendix  A  for  discussion  of  these  points).  Conse¬ 
quently  this  approach  results  in  a  larger  signal  to  noise  ratio, 
hereafter  denoted  by  S(u))/N(w),  for  all  frequencies. 

The  computational  procedure  used  for  obtaining  the  matched 
filter  spectral  estimates  as  well  as  the  conventional  estimates  is 
described  in  Appendix  C  and  shown  in  the  flow  chart  in  Figure  19. 
Briefly  the  procedure  is  as  follows: 

1.  The  filter  and  the  time  series  are  Fourier  transformed: 

y0(t)  -  Y0(u) 

x(t)  *  X(w) 

where  y  (t) *  the  filter,  is  the  large  reference  surface  wave  signal 
(obtained  empirically  or  generated  from  dispersion  curves  or  a 
combination  of  these). 

2.  The  Fourier  transform  of  the  time  series  is  multiplied  by 
the  complex  conjugate  of  the  reference  filter's  Fourier  transform  to 
give  the  Fourier  transform  of  the  matched  filtered  output: 


A(u>)  =  X (w) Y* (w) 

3.  This  product  is  transformed  back  to  the  time  domain: 

A(w)  -►  a ( t ) 
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4.  A  suitable  window  of  length  Tm  about  the  peak  in  a(t)  is 
selected  and  Fourier  transformed: 


a  ( t ) 
m 


A  (ui) 
m  ^  ‘ 


Note  that  A  (u>)  will  automatically  be  a  smoothed  spectral  estimate 
because  of  the  truncation  window  Tm  (Appendix  A) . 

S.  This  final  Fourier  transform  is  corrected  for  the  reference 
filter  response,  squared,  and  corrected  for  noise  just  as  in  the 
conventional  method 


IV“)I2 

|V0Cu)  i  - 


fVT„2)lS2<")|2 


(2) 


In  practice  it  is  both  convenient  and  advantageous  to  use  a  whitened 
reference  filter  such  that  i  YQ  ( to )  |  =  1  in  the  frequency  band  of 
interest  since  minima  in  Yq(u,)  were  found  to  produce  large  fluctua¬ 
tions  in  ps(<*>)  by  virtue  of  Am(w)  being  a  smoothed  spectral  estimate 
with  much  less  pronounced  minima  than  YQ(u).  Another  advantage  is 
that  the  matched  filtered  signal  should  be  of  shorter  duration  since 
it  is  a  band-limited  impulse  rather  than  an  autocorrelation;  this 
property  permits  use  of  a  shorter  signal  window.  That  is,  for  the 
whitened  case 


aw(t)  =  J  I  S  (w)  |  eiuitdu) 
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while  for  the  unwhitened  case 


a(t) 


(w)eiu,tda) 


■  r 


aw  (T) aW(t  +  T)dT • 


Thus  a(t)  is  approximately  twice  the  duration  of  aw(t)  since  it  is 
the  autocorrelation  of  aw(t). 

The  only  negative  feature  of  the  whitened  matched  filter  is 
that  as  a  detector  it  will  be  less  satisfactory  than  the  unwhitened 
filter.  However  we  are  not  addressing  the  detection  problem  in  this 
paper. 

Note  that  the  noise  correction  in  [2)  is  only  Tm/Tx  as  big  as 
that  in  equation  (1)  so  that  spectra  for  smaller  signals  should  be 
recoverable  at  lower  levels  than  for  the  direct  approach.  Also  by 
taking  a  finite  window  Tm  about  the  peak  we  are  automatically 
smoothing  the  signal  spectral  estimates.  The  noise  power  estimate 
in  (2)  will  be  the  same  as  that  in  (1)  if  the  noise  is  random  and 
stationary  (Papoulis,  1965  and  Appendix  A).  As  will  be  seen  later, 
this  assumption  about  the  noise  appears  to  be  warranted  for  the 
noise  samples  encountered  in  this  study. 

For  the  multichannel  (array)  spectral  estimates  three  different 
methods  were  employed  each  with  optional  weighting  among  channels. 
The  first  two  methods  (I  and  II)  consist  of  summing  with  weights 
the  individual  noise-corrected  spectra  obtained  from  equation  (1) 
and  from  equation  (2)  respectively.  This  gives  for  K  channels, 
suppressing  the  notation  for  frequency, 


K  K  _  K  ,  K  ,  K 

i?i  B‘  ■  tJi  B>  s‘  *  '  &  B‘ 


(3) 
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The  weights  are  typically  taken  to  be  inversely  proportional  to 
the  signal  to  noise  ratio  or  simply  inversely  proportional  to  the 
noise  RMS.  If  it  is  assumed  that  the  signal  spectrum  at  every 
station  is  the  same,  equation  (3)  becomes 


2  *i*3L/2  b. 
i-l  1  31  i-l  1 


SZ  ♦  I  B.  ni/  2  B. 
i-l  1  1  i*l  1 


s2  +  c 


Provided  the  noise  on  each  channel  is  random  and  stationary 
the  expected  value  of  the  term  £(w)  in  equations  (3)  and  (4)  is 
zero  since  for  each  we  overcorrect  for  noise  as  often  as  we 
undercorrect  in  equations  (1)  and  (2).  Moreover  the  variance  of 
C (w)  should  decrease  with  the  number  of  channels,  K.  Thus  the 
summing  will  give  an  improved  estimate  of  |S(u>)|2  with  a  variance 
in  the  case  of  uniform  weighting  given  by 


UHlalll  fi  .  Isk>I2i 

*  |N(u)  |  2  J 


(see  Appendix  B). 

The  third  technique  (Method  III)  is  to  beam  the  individual 
matched  filtered  time  series  and  compute  the  spectrum  of  the  matched 
filtered  signal  on  the  sum  trace  using  the  same  procedure  for  extract¬ 
ing  signal  spectra  as  used  for  the  individual  matched  filtered  traces 

(equation  2).  In  this  case  the  S/N  ratio  on  the  sum  trace  (time 

1/2 

domain)  is  approximately  (K  )  greater  than  on  individual  channels, 

2 

resulting  directly  in  a  better  spectral  estimate  of  S  (id)  with  a 
variance  in  the  case  of  uniform  weighting  given  by 


2  N(u] 


[i .  K  iiXsUii 
1  iN(«orJ 


(see  Appendix  B) . 
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Since  this  variance  is  always  smaller  than  for  the  other  two 
methods,  especially  at  small  signal  to  noise  ratios,  values  from 
Method  III  should  be  best  to  use.  As  before,  the  window  contain¬ 
ing  the  signal  remains  short  so  the  noise  contribution  is  minimum 
even  before  a  correction  for  it  is  made.  Also,  it  is  likely  that 
the  noise  on  the  sum  trace  wil]  be  more  neavly  stationary  than  on 
individual  channels  so  that  the  noise  correction  will  be  more  reli¬ 
able.  This  technique  has  the  further  advantage  that  only  one  spectral 
calculation  is  required  whereas  the  first  two  methods  each  require 
determination  of  K  individual  spectra.  This  latter  advantage  repre¬ 
sents  a  very  significant  saving  in  computer  time  and  makes  "on-line” 
spectral  calculations  feasible  to  consider  in  an  operational  system. 

These  methods  described  above  were  programmed  for  use  on  digital 
data.  The  computational  procedures  are  discussed  in  Appendix  C  and 
the  many  options  available  are  explained.  The  flow  chart  shown  in 
Figure  19  summarizes  the  program.  Basically  the  program  incorporates 
the  spectral  methods  (I,  II,  III)  into  a  previously  developed  pro¬ 
gram  (Alexander  and  Rabenstine,  1967)  for  matched  filtering  and 
array  processing  of  matched  filtered  seismograms. 

Note  that  if  the  signal  is  the  same  on  each  channel,  Method  III 
is  equivalent  to  beaming  with  a  constant  velocity,  and  then  match 
filtering  the  beam,  windowing,  and  transforming  to  get  the  spectra. 
This  would  further  reduce  the  calculations  for  Method  III  by  a 
factor  of  1/K.  The  size  of  error  incurred  at  LASA  by  this  approach 
is  unknown  and  should  be  investigated  before  a  real-time  system  is 
implemented. 

Each  of  the  above  techniques  should  yield  improved  estimates  of 
the  signal  spectrum  as  discussed  in  Appendix  B,  with  the  matched 
filter  beam  (Method  III)  preferred.  However  it  is  important  to 
investigate  which  approach  is  the  most  effective  in  practice.  We 
have  attempted  to  do  this  empirically  by  analyzing  test  cases  of 
actual  noise  where  the  signal  spectrum  is  known  exactly  and  by 
analyzing  actual  events  recorded  at  LASA.  The  results  of  these 
tests  are  presented  in  the  next  sections. 
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TEST  CASES 


Test  cases  were  generated  by  imbedding  a  set  of  synthetic 
Rayleigh  waves  into  real  noise  at  various  signal- to-noise  ratios. 

The  synthetic  signals  were  constructed  (Program  NEWSYN)  to 
correspond  to  surface  wave  signals  that  would  be  received  at  the 
nine  LASA  long  period  sites  (A,  E  and  F  rings  in  Figure  1A)  from 
an  epicenter  in  the  Greenland  Sea  (73. 4N,  6.8E).  With  the  instru¬ 
ment  amplitude  response  (Figure  IB)  included  all  had  identical 
power  spectra  as  shown  in  Figure  2A  and  phase  spectra  as  determined 
"rom  t lie  appropriate  epicentral  distances  (great  circle  path)  and 
the  dispersion  curve  shown  in  Figure  2B.  These  synthetic  seismograms 
are  shown  in  Figure  3, 

The  noise  was  real  LP  noise  recorded  digitally  at  5  samples 
per  second  at  the  LASA  on  March  30,  1967  and  reduced  to  1  sample 
per  second.  Typica1  normalized  power  spectra  of  these  noise  samples 
are  shown  in  Figure  4A. 

In  Appendix  A  it  is  shown  that  if  the  noise  is  random  and  sta¬ 
tionary,  matched  filtering  with  a  whitened  reference  filter  should 
not  change  the  noise  spectrum  on  the  matched  filtered  output.  (It 
also  follows  that  the  mean  and  RMS  in  the  time  domain  should  remain 
the  same.)  Figure  4B  shows  a  comparison  of  the  noise  power  spectra 
before  (solid  line)  and  after  (dashed  line)  matched  filtering  at 
LASA  stations.  A  30  minute  noise  window  was  used  throughout.  It  is 
clear  that  for  frequencies  lower  than  about  .05  Hz  the  power  spectra 
are  essentially  unchanged  in  both  level  and  shape.  The  shape  remains 
similar  at  the  shorter  periods  but  the  level  may  differ  by  up  to  a 
factor  of  2.  The  results  shown  in  Figure  4B  are  typical  of  the  other 
LASA  channels  and  other  noise  cases  tested.  We  conclude,  therefore, 
that  at  LASA  the  noise  on  each  channel  is  sufficiently  random  and 
stationary  for  frequencies  below  .05  Hz  that  matched  filtering  does 
not  appreciably  change  the  noise  spectrum.  However  for  higher 
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frequencies  the  stationarity  assumption  appears  to  bo  only  approxi¬ 
mately  valid. 

Also  shown  in  Figure  4B  is  the  sum  trace  noise  spectrum  for 
the  3  channels  shown;  its  power  spectral  level  is  approximately 
3  times  smaller  than  the  individual  channel  levels.  This  result 
implies  that  most  of  the  LP  noise  is  not  correlated  from  station 
to  station  at  LASA.  Alexander  and  Rabenstine,  (1967),  reached  a 
similar  conclusion,  and  all  the  results  in  this  study  using  more 
stations  give  noise  reductions  on  summing  appropriate  to  uncorre¬ 
lated  noise,  provided  the  stations  are  20-2S  km  or  more  apart. 

Figure  4C  shows  the  effect  of  window  length  on  the  matched 
filter  spectral  estimates  for  the  signal  spectrum  in  the  synthetic 
test  cases.  The  window  lengths  are  shown  at  the  top  centered  on  the 
matched  filtered  signal  (whitened  reference)  which  has  a  true  norma¬ 
lized  spectrum  given  by  the  dashed  curve  in  the  power  spectrum  plots 
in  the  same  Figure.  It  is  evident  that  a  window  of  300  seconds  length 
or  more  does  not  significantly  distort  the  signal  spectrum  in  this 
case,  whereas  the  50  and  100  second  windows  result  in  over- smoothing . 
Thus  the  shortest  window  consistent  witli  suitably  undistorted  signal 
spectra  is  300  seconds;  this  window  length  was  used  throughout  the 
analysis  of  test  cases. 

Test  case  seismograms 

The  signal  to  noise  ratio  used  in  preparing  the  test  cases  were 
defined  as  1/2  the  peak  to  peak  value  of  the  detrended  signal  to  *he 
RMS  amplitude  of  the  detrended  noise.  This  gave  a  value  of  input 
signal  to  noise  which  was  probably  conservative  (too  high)  by  3  db 
since  the  1/2  peak  to  peak  value  is  typically  about  1.4  times  the 
RMS  amplitude  of  the  signal  (assuming  a  long  dispersed  signal  wave¬ 
form).  However,  this  definition  of  S/N  allows  comparison  with 
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previous  studies  (Alexander  and  Rahenstine,  SDL  Report  Nos.  17S 
and  194,  1967)  where  a  similar  definition  has  been  used.  Figure 
5A  shows  a  synthetic  case  where  the  signal  is  mixed  with  the  noise 
at  S/N  »  8.  The  signal,  the  noise,  signal  ♦  noise,  and  the  matched 
filtered  traces  are  plotted  with  a  common  vertical  scale  to  illuus- 
trate  graphically  the  meaning  of  S/N  ■  8  as  defined  above  and  to 
show  that  the  matched  filtering  does  not  change  the  RMS  noise 
level.  The  whitened  matched  filter  option  was  used  in  this  case 
and  throughout  the  remainder  of  this  study,  because  as  discussed 
earlier  it  generally  produces  signals  of  shorter  duration  than 
the  autocorrelation  resulting  from  standard  matched  filtering 
(i.e.  direct  cross-correlation  of  reference  signal  with  the  seismo¬ 
gram),  and  because  we  found  that  minima  in  the  reference  signal 
spectrum  caused  fluctuations  in  the  estimate  of  the  weak  signal 
spectrum.  These  spurious  fluctuations  are  due  both  to  noise 
"filling  up"  any  actual  holes  in  the  weak  signal  spectrum  and  to 
the  inherent  smoothing  from  using  a  short  signal  window  on  the 
matched  filter  output  (Figure  4C). 

Figures  5B  and  SC  show  the  individual  input  traces  for  S/N-8 
and  the  corresponding  matched  filtered  outputs.  The  measured  S/N 
on  each  matched  filtered  trace  is  indicated.  Figure  5C  also  shows 
the  phased  sum  (beam)  of  nine  individual  matched  filtered  outputs. 
Note  that  the  matched  filtering  results  in  individual  S/N  ratio 
gains  in  the  time  domain  of  slightly  over  two  while  the  sum  of 
nine  elements  results  in  a  factor  of  over  seven  S/N  ratio  gain 
above  the  original  individual  inputs.  These  results  are  consistent 
with  the  earlier  results  of  Alexander  and  Rabenstine,  (1967)  for 
LASA. 

Figure  6A  shows  the  test  case  synthesis  for  S/N  *  4  again  with 
all  traces  plotted  with  a  common  amplitude  scale.  Figures  6B  and 
6C  show  the  individual  input  traces  and  matched  filtered  outputs. 
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Figure  6C  shows  the  sum  of  nine  matched  filtered  outputs.  As 
expected  we  achieve  somewhat  over  a  factor  of  2  S/N  gain  on 
matched  filtered  outputs  and  about  a  factor  of  seven  on  the  sum. 

Figure  7A  shows  the  test  case  synthesis  for  S/N  -  2.  Here  the 
original  signal  is  not  evident  in  the  noise  whereas  it  is  in  the 
matched  filtered  output.  Figures  7B  and  7C  show  the  input  traces 
and  matched  filtered  traces  for  this  case.  In  most  instances  the 
signal  is  evident  on  the  individual  matched  filtered  traces  but  is 
nearly  at  the  threshold  where  many  false  alarms  would  occur  in  a 
search  for  signals  at  this  level.  The  array  sum  of  nine  elements, 
however,  has  an  unmistakeable  signal  present  (bottom  trace  on  7C) . 

Figure  8A  shows  the  test  case  synthesis  for  S/N  ■  1.  In  this 
case  the  signal  is  not  evident  on  the  signal-to-noise  trace  nor  is 
it  evident  on  the  matched  filtered  output.  Thus,  the  threshhold  for 
single  channel  detection  (assuming  identical  signals)  is  between 
that  shown  in  Figure  7A  and  that  shown  in  8A.  Figures  8B  and  8C 
verify  this  conclusion  but  show  that  the  threshold  for  detection 
is  below  S/N  -  1  for  nine  LASA  channels  since  the  signal  on  the  sum 
trace  in  Figure  8C  is  unmistakeable.  We  further  verified  this 
empirical  result  by  imbedding  the  signal  at  various  positions  and 
repeating  the  analysis. 

In  principle  one  expects  individual  (time  domain)  gains  in  S/N 
resulting  from  matched  filtering  to  be  proportional  to  the  time- 
bandwidth  product  of  the  original  signal  and  the  gain  on  summing  to 
be  /TT  where  K  is  the  number  of  sensors,  so  the  reader  is  cautioned 
to  generalize  the  above  results  to  other  signals  and  numbers  of 
stations  with  these  constraints  in  mind. 

Test  case  spectra 

We  now  turn  to  a  discussion  of  the  various  spectral  estimates 
obtained  from  analyzing  these  test  cases.  By  comparing  typi  ;al  single 
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station  spectral  estimates  and  their  sums  with  the  test  seismograms 
discussed  above  one  can  begin  to  get  a  feel  for  what  spectral  esti* 
mates  to  believe  vis  &  vis  the  actual  seismograms  from  which  they 
were  derived. 

Figure  9  shows  spectra  for  five  of  nine  individual  LASA  stations 
as  indicated  for  the  S/N  ■  8  case  (see  corresponding  seismograms  in 
Figures  SB  and  SC).  The  dashed  black  line  is  the  input  spectrum, 
the  black  line  is  the  noise* corrected  spectrum  obtained  by  conven¬ 
tional  analysis  (equation  (1)),  and  the  red  line  is  the  noise- 
corrected  matched  filter  signal  spectrum  (equation  (2)).  For  the  sum 
the  dashed  line  shows  the  actual  spectrum,  the  black  line  the  sum  of 
nine  conventionally  computed  spectra  (equation  (3)),  the  solid  red 
line  the  sum  of  the  nine  individual  matched  filter  spectra,  and  the 
dashed  red  line  the  noise-corrected  spectrum  from  the  matched  filtered 
sum  trace  shown  in  Figure  SC.  The  corresponding  noise  power  spectra 
for  these  individual  channels  are  shown  in  Figure  4A.  It  is  clear 
that  the  conventional  and  matched  filter  estimates  are  very  similar 
to  one  another  and  that  the  sums  all  give  much  better  representation 
of  the  actual  spectrum  than  do  individual  ones  even  for  this  rela¬ 
tively  large  S/N  case.  As  indicated  earlier  the  matched  filter  spec¬ 
trum  should  be  a  smoothed  version  of  the  actual  signal  spectrum  and 
this  effect  can  be  seen  by  comparing  the  solid  red  and  black  lines 
in  Figure  9.  The  signal  window  used  was  300  seconds.  The  conclusion 
is  that  at  this  S/N  level  the  methods  are  all  comparable  in  perfor¬ 
mance  and  give  good  spectral  estimates,  especially  the  sums. 

Figure  10  shows  the  same  display  as  does  Figure  9  but  for  S/N*4 
(see  Figures  6B  and  6C  for  corresponding  seismograms).  In  this  case 
the  individual  matched  filter  spectra  (solid  red  lines)  are  as  good 
or  better  than  the  conventional  ones  (solid  black  lines)  and  the 
simmed  spectra  are  better  than  any  individual  spectral  estimate.  Again, 
the  spectrum  from  the  matched  filter  sum  (Figure  6C)  is  the  best  in 
this  case  as  is  expected  theoretically  (Appendix  B) .  Thus  at  about 
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S/N  ■  4  the  matched  filter  spectra  begin  to  give  more  reliable 
spectral  estimates  than  do  the  conventional  ones,  although  the 
conventional  sum  is  still  quite  acceptable  in  level  and  shape. 

Figure  11  shows  the  same  display  as  Figures  9  and  10  hut 
for  S/N  ■  2  (see  Figures  ?B  and  7C  for  corresponding  seismograms ) . 

In  this  case  the  individual  spectra  are  all  significantly  in  error, 
although  the  matched  filter  spectra  define  the  location  of  the 
actual  spectral  peak  somewhat  more  reliably  than  do  the  conventional 
estimates.  However,  sums  of  the  individual  matched  filter  spectra 
and  the  spectrum  from  the  matched  filter  sum  trace  (Figure  7C)  give 
reliable  spectral  values  and  shape  around  the  peak  and  give  reason¬ 
ably  accurate  signal  levels  at  the  higher  frequencies.  The  sum  of 
conventional  spectra  (solid  black  line)  is  significantly  in  error. 
Again,  the  best  result  comes  from  the  matched  filter  sum  trace. 

Finally  Figure  12  shows  the  results  for  S/N  •  1  (see  Figures  8B 
and  8C  for  the  corresponding  seismograms).  In  this  case  none  of  the 
individual  spectra  comes  close  to  representing  the  actual  spectrum. 
Only  the  sum  of  individual  matched  filter  spectra  (solid  red  line) 
and  the  spectrum  of  the  matched  filter  sum  trace  (Figure  8C)  produce 
reasonable  spectral  estimates,  and  even  then  the  spectral  shape  is 
significantly  distorted,  except  at  the  peak. 

The  performance  of  the  various  methods  can  be  assessed  by 
comparing  quantitatively  the  array  estimates  obtained  in  the  test 
cases  with  the  actual  signal  power  in  the  array  sums.  Figure  13A 
shows  the  variance  i.e.  (estimate  -  actual)2  as  a  function  of 
frequency  for  S/N  -  4  using  nine  LASA  elements  (see  Figure  10  for 
the  corresponding  spectral  estimates).  Clearly  the  most  accurate 
result  (smallest  variance)  overall  is  given  by  Method  III  and  the 
least  accurate  by  Method  I,  with  Method  II  between  them  but  still 
much  better  than  Method  I.  Figure  13B  shows  the  variance  with 
frequency  for  S/N  ■  Z.  (See  Figure  11  for  the  corresponding  spectral 
estimates.)  Finally  Figure  13C  shows  the  variance  with  frequency 
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for  S/N  =  1.  (See  Figure  12  for  the  corresponding  spectral  esti¬ 
mates.)  Method  III  is  consistently  the  most  accurate  in  each 
instance.  However,  as  expected,  Method  II  is  also  much  better 
than  Method  I  because  the  noise  power  is  smaller  by  the  ratio 
of  the  respective  window  lengths  taken.  These  results  are  in 
accord  with  what  is  expected  from  the  simplified  theory  given 
in  Appendix  B.  There  it  is  shown  that  Method  III  is  the  most 
accurate  and  should  get  progressively  better  than  the  other 
sums  as  the  S/N  decreases. 

Since  the  slowly  varying  signal  power  spectrum  was  norma¬ 
lized  to  its  maximum  and  kept  fixed  for  all  S/N  cases  the  square 
root  of  the  variance  values  times  100  represents  approximately 
the  percentage  error  in  the  spectral  estimates  in  each  case. 

Thus  from  Figures  13A,  B,  C  it  is  clear  that  Method  I  breaks  down 
(i.e.  errors  of  over  100  percent  at  some  frequencies)  between 
S/N  =  4  and  S/N  =  2;  Method  II  breaks  down  I etween  S/N  =  2  and 
S/N  =  1;  on  the  other  hand  Method  III  at  S/N  =  1  has  a  maximum 
error  for  the  entire  frequency  band  of  about  70  percent  and  an 
error  at  the  peak  of  the  signal  of  only  10  percent  (variance  of 
1  percent).  It  must  be  cautioned  that  these  numerical  results 
apply  only  to  the  particular  test  case  signals  and  noise  samples 
used  with  sums  over  nine  LASA  elements.  Therefore,  they  cannot 
be  taken  as  representative  of  the  errors  expected  in  general  from 
applying  these  techniques,  although  by  design  the  test  cases 
should  be  realistic  for  signals  from  the  Greenland  Sea  area 
observed  at  LASa. 

Thus  the  practical  threshhold  for  the  spectral  methods  con¬ 
sidered  here  is  between  S/N  of  1  and  2  for  nine  LASA  channels. 
The  best  estimator  both  theoretically  and  observationally  seems 
to  be  the  noise-corrected  spectrum  of  the  matched  filter  sum 
tracc  signal  (Method  III).  Certainly  it  is  consistently  at  lease 
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as  good  as  either  of  the  other  sums.  This  is  an  important  finding, 
because  arrays  such  as  LASA  or  even  larger  continental  sized  arrays 
can  be  used  to  obtain  reliable  spectra  for  weak  events  without 
having  to  compute  the  individual  station  spectra. 

Again  it  must  be  stressed  that  this  particular  threshhold  S/N 
applies  for  LASA  (nine  elements)  for  a  particular  signal  duration 
and  bandwidth.  While  this  case  may  be  roughly  representative  of 
teleseismic  Rayleigh  waves  received  at  LASA  an  empirical  threshhold 
test  similar  to  that  described  above  should  be  carried  out  for  ref¬ 
erence  events  from  each  source  region  of  interest.  However,  we 
still  can  conclude  that  this  experiment  has  demonstrated  that  reli¬ 
able  spectra  for  realistic  weak  surface  wave  signals  can  be  obtained 
using  the  array  methods  discussed  above,  and  that  Method  III  is  prefer 
erable  to  use  from  the  standpoint  of  both  accuracy  and  economy  in 


APPLICATION  TO  ACTUAL  EVENTS 


Using  the  methods  presented  earlier  and  the  results  of  the 
test  cases,  actual  events  were  analyzed  and  the  three  spectral 
estimates  compared.  We  will  discuss  only  two  events  in  detail 
here  because  a  separate  report  will  be  published  presenting  many 
observed  spectra  for  various  regions. 

Epicenter  data  for  the  events  used  are  given  in  Table  1. 

For  the  Greenland  Sea  region  Figures  14A,  14B,  and  14C  show 
the  reference  event,  the  input  trace,  and  the  matched  filtered 
output  respectively  for  each  of  nine  LASA  stations  (Z  component)  in 
the  D,  E,  and  F  rings.  Figure  14C  also  shows  the  sum  trace.  Each 
trace  is  normalized  to  its  maximum  value  in  Figures  14A-C.  Because 
the  noise  RMS  is  unchanged  by  matched  filtering,  the  apparent  reduc¬ 
tion  in  noise  level  shown  for  individual  stations  in  this  figure  is 
a  measure  of  the  minimum  S/N  gain  achieved  in  each  case.  The  matched 
filtered  traces'  S/N  levels  are  between  those  for  the  test  cases  at 
S/N  =  4  and  S/N  *  2  (Figures  6C  and  7C  respectively)  or  approxi¬ 
mately  at  S/N  *  3.  Thus  the  spectra  obtained  by  each  of  the  three 
methods  should  be  comparable  in  accuracy  to  those  for  the  test  cases 
at  S/N  =  4  and  S/N  =  2  (Figures  10  and  11). 

Figures  ISA,  15B,  and  15C  show  spectral  results  for  three  of  the 
nine  individual  LASA  channels.  Figures  14A,  14B,  and  14C  show  the 
corresponding  time  series.  In  each  case  we  show  the  reference  event 
spectrum,  the  conventional  S  +  N  spectrum,  the  noise  spectrum,  the 
MF  spectrum  (S  ♦  N) ,  and  the  noise-corrected  conventional  and  MF 
spectra.  Since  the  instrument  responses  are  nominally  the  same  for 
all  the  LASA  sites  (Figure  IB)  and  unchanged  from  event  to  event  we 
do  net  need  to  make  instrument  corrections  to  obtain  valid  comparisons 
of  the  spectral  estimates.  Both  the  conventional  and  matched  filtered 
results  at  the  single  stations  give  very  similar  spectra  in  both 
level  and  shape,  with  the  matched  filtered  spectra  smoother,  as 
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expected.  Note  the  similarity  in  spectral  shape  between  the  refer* 
ence  event  spectrum  and  the  smaller  event  spectrum  even  though 
the  signal  power  everywhere  differs  by  a  factor  of  over  100  (and 
surface  wave  magnitude  by  more  than  1  unit  as  compared  with  the 
USC§GS  body-wave  magnitudes  (Table  I)  which  differ  by  only  .3  units, at 
4.9  vs  4.6.  This  remarkable  similarity  in  spectral  shape  implies 
that  the  source  mechanism  for  both  events  was  essentially  the  same. 
This  particular  pair  of  events  also  illustrates  the  importance  of 
using  a  rather  broad  frequency  band  for  estimating  surface  wave 
magnitude.  From  Figures  15A,  1SB,  and  15C  the  signal  power  at  20 
seconds  period  where  conventional  surface  wave  magnitudes  are  deter¬ 
mined  is  almost  an  order  of  magnitude  smaller  than  the  peak  power  at 
around  30  seconds  period.  Also  the  noise  contamination  is  greatest 
at  around  20  seconds  period  (Figures  15A,  B  and  C) . 

Figure  15D  compares  the  different  array  summed  spectra  for  the 
A,  E,  and  F  rings  (nine  elements)  at  l.ASA  for  the  small  Greenland  Sea 
event  just  discussed.  It  is  clear  that  all  three  methods  give 
essentially  the  same  signal  power  estimates.  By  analogy  with  the 
corresponding  S/N  test  case  (Figures  9  and  10)  we  conclude  that  these 
are  reliable  spectral  determinations.  Also  the  consistency  among  the 
three  differer  estimates  of  the  spectrum  suggest  that  it  is  a 

reliable  rest 

It  is  imp  to  note  that  for  this  real  case  the  previous 

conclusion  that  Ned  od  III  is  the  best  to  use  is  supported,  because 
the  single  estimate  from  the  matched  filter  beam  trace  is  very 
nearly  the  same  as  that  produced  by  the  other  two  array  methods,  and 
requires  less  computing  time. 

For  the  Yunnan  events  in  Table  I  we  show  in  Figures  16A,  16B, 
and  16C  the  reference  event  (  23  September  1966)  the  input  trace  for 
29  March  1967  and  the  matched  filtered  output  respectively  for  each 
of  nine  LASA  stations  in  the  AO,  0,  !:,  and  F  rings.  Figure  16C  also 
shows  the  matched  filtered  sum  trace.  Each  trace  is  normalized  to  its 
maximum  value  in  Figures  16A,  B  and  C.  Because  the  noise  RMS  is 
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unchanged  by  matched  filtering  the  apparent  reduction  in  noise 
level  in  each  case  is  a  measure  of  the  minimum  S/N  gain  achieved 
just  as  in  the  Greenland  Sea  event  discussed  above.  The  matched 
filtered  traces  S/N  levels  are  between  those  for  the  test  case 
at  S/N  ■  4  and  S/N  *  2  (Figures  6C  and  7C  respectively)  but  some¬ 
what  closer  to  the  S/N  *  4  case.  Thus  the  spectra  obtained  by  each 
of  the  three  methods  should  be  comparable  in  accuracy  to  those  for 
the  test  cases  at  S/N  *  4  and  S/N  *  2  (Figures  10  and  11)  but 
closer  to  the  case  in  Figure  10. 

Figures  17A,  17B,  and  17C  show  the  spectral  results  for  three 
of  the  LASA  stations  used  for  the  Yunnan  event  (refer  to  Figures 
16A,  16B,  and  16C  for  the  corresponding  time  series).  These  are 
typical  of  the  single-station  results  for  this  event.  As  before 
we  show  in  each  case  the  reference  event  spectrum,  the  conventional 
S  +  N  spectrum,  the  noise  spectrum,  the  MF  spectrum  S  +  N,  and 
the  noise-corrected  conventional  and  MF  spectra.  Again  the  conven¬ 
tional  and  matched  filtered  results  at  the  single  stations  are  very 
similar  in  both  level  and  shape ,  wi  th  the  matched  filter  spectra 
somewhat  smoother  as  expected  from  theory  and  the  test  cases  at 
S/N  =  4.  In  contrast  to  the  Greenland  Sea  events  discussed  above 
the  small  event  and  the  reference  event  spectra  match  only  approxi¬ 
mately  in  shape  with  the  smaller  event  enriched  in  higher  frequencies 
compared  to  the  reference.  Also  there  appears  to  be  somewhat  more 
variability  in  shape  and  level  from  station  to  station  over  LASA  for 
the  Yunnan  pair  as  compared  to  the  Greenland  Sea  pair  (Figures 
15A,  B,  and  C) ;  however  the  shapes  and  levels  are  still  not  greatly 
different  from  station  to  station. 

Figures  17D  compares  the  different  array  summed  spectra  using 
Methods  I,  II,  and  III  on  channels  at  LASA  for  the  small  Yunnan  event 
(Table  I).  It  is  clear  that  all  three  methods  give  essentially  the 
same  signal  power  estimates.  By  analogy  with  the  corresponding  S/N 
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test  case  (Figure  9)  we  conclude  that  these  are  reliable  spectral 
determinations  comparable  in  accuracy  or  better  than  the  estimates 
for  the  small  Greenland  sea  event  of  18  November  1966  (Table  1  and 
Figure  1SD)  discussed  earlier.  Again  the  consistency  among  all  the 
estimates  further  suggests  that  the  spectra  are  reliable.  Also  we 
emphasize  that  the  single  estimate  from  the  matched  filtered  beam 
trace  (Method  III)  is  in  excellent  agreement  with  the  other  spectral 
estimates  for  this  real  case.  This  is  to  be  expected  theoretically 
and  by  analogy  with  results  for  the  test  cases  at  S/N  s  4.  Therefore 
our  conclusion  that  Method  III  is  preferable  to  use  operationally 
because  of  the  computational  savings  seems  to  hold  in  practice. 

While  the  two  actual  events  discussed  hardly  allows  us  to 
generalize  we  have  demonstrated  that  reliable  spectra  can  readily 
be  obtained  teleseismically  using  any  of  the  three  methods  (I,  II, 
or  III)  for  events  of  m^  approaching  4.5.  By  comparison  with  the 
test  cases  we  can  infer  that  Method  III  should  give  reasonably 
good  spectral  estimates  at  S/N  ratios  approaching  1  (assuming  a 
nine  element  array)  at  LASA  so  a  reasonable  expected  threshhold 
would  be  well  below  m^  *  4.S  if  the  examples  used  here  are  at  all 
typical.  Use  of  more  than  nine  elements  in  the  array  will  lower 
the  threshhold  still  more.  This  question  of  threshhold  for  spectral 
determinations  will  be  deferred  to  a  following  report  now  in  prepara¬ 
tion  summarizing  spectral  results  for  many  different  events  and 
different  source  regions. 
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CONCLUSIONS 


Based  on  the  findings  of  this  investigation  we  conclude  that: 

1.  Reliable  signal  spectra  can  be  obtained  down  to  values  of 
S/N  ratio  approaching  unity  by  combined  use  of  matched  filtering 
and  array  processing  of  9  LASA  channels.  In  principle  the  threshold 
can  be  further  lowered  simply  by  using  more  sensors. 

2.  The  "best"  spectral  estimator  both  theoretically  and  experi¬ 
mentally  consists  of  first  beaming  individual  matched  filtered  seis¬ 
mograms  and  then  computing  a  single  noise- corrected  spectrum  from 
this  beamed  trace.  This  result  is  extremely  important  from  a  practical 
standpoint  since  signficant  computing  time  is  saved  by  having  to 
compute  only  one  noise-corrected  spectrum. 

3.  It  is  feasible  to  consider  "on  line"  spectral  calculations 
since  only  the  matched  filtered  beams  are  required  using  the  techniques 
developed  here. 

4.  A  signal  window  of  about  300  seconds  on  the  matched  filtered 
trace  was  shown  to  be  a  good  practical  length  to  use  so  as  to  pre¬ 
serve  spectral  detail  of  signals  in  the  15-S0  second  period  range. 

At  the  same  time  the  300  second  window  reduces  the  total  noise  contami¬ 
nation  compared  to  signal  windows  required  for  long,  dispersed  wave 
trains . 

5.  Using  the  techniques  and  computer  programs  developed  in  this 
study  surface  wave  spectra  can  he  routinely  obtained  and  saved  for 
many  events  from  source  regions  of  interest. 

6.  It  is  both  desirable  and  convenient  to  use  whitened  refer¬ 
ence  signals  because  (a)  the  resulting  time  domain  signals  (band 
limited  impulses)  are  of  shorter  duration  than  cross  correlations 
resulting  from  standard  matched  filtering  thus  allowing  use  of  shorter 
signal  windows  and  (b)  spurious  peaks  in  the  weak  signal  spectral 
estimates  resulting  from  minima  in  the  reference  signal  spectrum  are 
avoided . 
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Figure  1A.  Configuration  of  the  Large  Aperture  Seismic  Array 
in  Montana. 


RELATIVE  MAGNIF1  CATION 


Figure  IB.  LASA  long  period  system  "A"  response 


Figure  2A.  Normalized  power  spectrum  of  synthetic  signals  computed 
for  the  LASA,  A,  t,  and  F  rings. 
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Figure  4A.  Typical  normalized  noise  power  spectra  at  six  LASA 

stations  for  a  30  minute  noise  sample  on  30  March  1967 
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Figure  4B.  Comparison  of  directly  computed  noise  power  spectra 
and  noise  power  spectra  from  the  matched  filtered 
output  at  three  LASA  stations. 
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Figure  4C.  Effect  of  window  length  on  matched  filtered  signal 
spectral  estimates. 


Figure  5A.  Typical  synthetic  test  case  normalized  to  noise  RMS 
showing  to  scale  the  mixing  of  signal  and  noise  and 
the  matched  filtered  output  for  S/N  =  8. 


Figure  5B.  Typical  synthetic  test  case  showing  input  traces 
matched  filtered  outputs  normalized  to  peak  ampli 
for  four  LASA  stations.  S/N  =  8. 
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Figure  6C.  Typical  synthetic  test  case  showing  input  traces  a 
matched  filtered  outputs  normalized  to  peak  amplit 
for  four  LASA  stations  and  the  phased  sum  trace  of 
nine  matched  filtered  outputs.  S/N  =  4. 
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Figure  7A.  Typical  synthetic  test  case  normalized  to  noise  RMS 
showing  to  scale  the  mixing  of  signal  and  noise  and 
the  matched  filtered  output.  Amplitude  scales  are 
the  same  on  all  traces.  S/N  =  2. 


matched  filtered  outputs  normalized  to  peak  amplitude 
for  four  LASA  stations.  S/N  =  2. 
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Figure  8A.  Typical  synthetic  test  case  normalized  to  noise  RMS 
showing  to  scale  the  mixing  of  signal  and  noise  and 
the  matched  filtered  output.  Amplitude  scales  are  the 
same  on  all  traces.  S/N  =  1. 
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Figure  9.  Synthetic  test  case  spectra  corrected  for  noise  for 
five  LASA  stations  and  the  sums  of  nine  LASA  station 
spectra  compared  with  the  original  signal  spectrum. 
S/N  =  8. 
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Figure  10.  Synthetic  test  case  spectra  corrected  for  noise  for 
five  LASA  stations  and  the  sums  of  nine  LASA  station 
spectra  compared  with  the  original  signal  spectrum. 
S/N  =  4. 
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Figure  11.  Synthetic  test  case  spectra  corrected  for  noise  for 
five  LASA  stations  and  the  sums  of  nine  LASA  station 
spectra  compared  with  the  original  signal  spectrum. 
S/N  =  2. 
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Figure  12.  Synthetic  test  case  spectra  corrected  for  noise  for 
five  LASA  stations  and  the  sums  of  nine  LASA  station 
spectra  compared  with  the  original  signal  spectrum. 
S/N  =  1. 
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Figure  ISA.  Individual  LASA  station  spectra  for  Greenland  Sea 

events,  trace  D4Z  uncorrected  for  instrument  respons 
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Figure  15B,  Individual  LASA  station  spectra  for  Greenland  Sea 

events,  trace  H4Z  uncorrected  for  instrument  response 


Figure  ISC.  Individual  LASA  station  spectra  for  Greenland  Sea 

events,  trace  I- 3 Z  uncorrccted  for  instrument  response 
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Figure  151).  Summed  spectra  and  sum  trace  spectra  of  nine  LASA 
stations  for  Greenland  Sea  events,  uncorrectcd  for 
instrument  response. 
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Figure  17A. 


Individual  LASA  station  spectra  for  Yunnan 
events,  trace  F4Z  uncorrected  for  instrument 
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Figure  17B. 


Individual  LASA  station  spectra  for  Yunnan,  China 
events,  trace  E3Z  uncorrected  for  instrument  response. 


Figure  17C.  Individual  LASA  station  spectra  for  Yunnan,  China 

events,  trace  E2Z  uncorrected  for  instrument  response 
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APPENDIX  A 

MATCHED  FILTER  METHOD  AND  ITS 
ASSOCIATED  SIGNAL  TO  NOISE  GAIN 


We  present  here  a  discussion  of  the  details  of  the  matched 
filter  operation  which  shows  that  because  the  signal  is  compressed 
in  time  while  the  noise  is  not,  S/N  is  increased.  Basically  the 
presentation  is  the  same  as  that  in  the  Appendix  [Alexander  and 
Rabenstine,  1967)  with  some  additions  and  revisions  relevant  to 
this  study. 

The  test  seismogram  x(t)  has  the  Fourier  representation 
x  (t )  -  J  X  (u)  eiu)tdw  (A- 1) 

-  OO 

and  energy  (by  Parseval's  identity) 

Ex  *  J  |  X  (ui)  |  2  du  (A- 2) 

The  matched  filter  has  the  Fourier  representation 

100 

X(w)Y*(w)  elwtda)  (A-3) 

-  00 

where  Y£(u)  is  the  complex  conjugate  of  the  Fourier  transform  of 
the  reference  signal  y(t),  (Y  (w)  «  | Yq (u) | e*9 ^ ) .  The  associated 
energy  of  the  matched  filter  output  is 

E  *  {  |X(w)  |  2  I  Yq (u>)  |  2  du  (A-4) 

•  00 

This  is  equivalent  to  the  energy  in  x(t)  after  it  is  filtered 

with  a  function  whose  spectrum  is  Y  (u) .  Note  that  if  we  whiten 

Y  (w)  (i.e.  divide  by  jY  (w)|)  so  that  instead  of  Y*(u)  we  have 
o  - i a  rwi  ^  ® 

Y*/|Y  (u)|  *  e  1  J  in  equation  (A-3),  this  equation  becomes 

aw(t)  *  J  X(oo)e"*£9  (uj-utjd^  (A-5) 

-  <X> 


-Al- 


with  energy 


En  “1  I X(“) I 2dw  (A- 6) 

Thus,  if  the  reference  signal  is  whitened  before  applying  it 

as  a  matched  filter  the  energy  of  x(t)  is  conserved  in  the  operation 

since  E  =  E  . 
aw  x 

What  is  left  to  show  is  that  the  signal  is  compressed  in  time 
while  the  noise  is  not. 

If  we  take  x(t)  =  s(t)  +  n(t)  (signal  +  random  noise),  then, 


X(w)  =  |  S  (u>)  |  exp  £i($os(u>)  -  +  (w))j 

+  |  N(w)  |  exp  [i(V(oi)  ♦  4>i(tD))J 


(A- 7) 


Vu) 


and 


where 


|Y0(-)|  exp  [-i(*oy(M)  - 
|Y0(w)t  e'i9(u)) 


(A-  8 ) 


4>  (u) 

os  v 

<t>  (u) 

(oj) 

C  (co) 


A 

|S(u)  | 


phase  spectrum  at  the  source  for  s(t) 
phase  spectrum  at  the  source  for  y(t) 

instrument  phase  spectrum  (assumed  identical  for  x  and  y) 
phase  velocity  spectrum 
epicentral  distance 
amplitude  spectrum  of  s(t) 


-A2- 


I Y  (  u>)  |  -  amplitude  spectrum  of  y(t) 

0 (u>)  =  total  phase  associated  with  y(t) 

't'(w)  =  random  phase  associated  with  n(t) 

therefore , 

XY*  =  |  S  (u»)  |  |Y0(u)|  exp[i(^0S(W)  -  *oy(w))] 
♦  !  N  (to)  |  |YQ(w)|  exp[i(Y(w)  -  6(w))j 


(A- 9) 


The  first  term  on  the  right  hand  side  of  the  equation  (A -  9 )  is 
the  matched  filter  signal  spectrum  and  the  second  term  is  the  noise 
spectrum.  We  will  discuss  each  term  separately.  From  the  first  term, 
the  signal  on  the  matched  filter  seismogram  is  given  by 

500 

|S(W)|  |Y0(a))j  exp[i  (4>os(w)  -  <t>oy(w)  ♦  «t)]  dw  (A- 10 ) 


If  both  sit)  and  y(t)  have  the  same  phase  spectrum  at  the  source 


then  <t>os  *  4>oy  and  (A- 10)  reduces  to 


sf(t)  =  J  |S(u»)  |  ! Y  (oo)  i  e1 


iwt 


doo 


(A- 11) 


or  for  the  whitened  reference  signal  case  (i.e.  | YQ (to)  j  =  1) 


’fw 


(t)  =  j  I  S  (co)  |  eiwt 


dw 


(A- 12) 


This  is  the  band-limited  impulse  representing  the  signal  of  shortest 
duration  which  has  the  spectrum  of  the  signal.  Therefore  this  signal 
also  lias  the  largest  possible  energy  density  in  the  time  domain.  To 
demonstrate  this  we  compare  the  energy  density  of  s^w(t)  with  that 


-A3- 


of  the  dispersed  signal  as  recorded 

s(t)  =  J  |S(io)|  exp[i  (0OS  (w)  -  •jj—j  *  +  wt)]  du)  (A-13) 

For  s^w(t)  the  maximum  energy  density  occurs  at  t  =  0,  since  it 
is  here  that  all  the  frequencies  add  together  in  phase  (zero  group 
delay  for  all  frequencies)  and 

I£»  «  00 

|  S  ( oj )  ]  dto  j  |  S  (id)  |  elaJtdu)  for  all  t  /  0, 

-  00  -  00 

Note  that  S£w(t)  is  not  "causal"  because  it  is  non-zero  earlier  than 
t  =  0.  That  is,  s  £- w  ( t )  is  symmetric  about  t  =  0  and  represents  what 
is  obtained  by  band-pass  filtering  a  delta  function  with  a  phaseless 
filter  of  spectrum  |S(w)|.  To  make  it  causal  we  would  construct  the 
minimum  phase  wavelet  with  the  spectrum  |  S  (to)  | .  For  s(t)  the  energy 
is  spread  over  a  longer  time  window  because  the  group  delays  (energy 
arrival  times)  are  frequency-dependent;  that  is,  the  energy  arrival 
at  frequency  uj  occurs  at  a  value  of  t  such  that 


d 

ITS 


+  *i(u)  +  wt] =  *os(u,)  *  ufcr  +  1  *  0 


or 


t 


£ 

uT^T 


4>qS(w)  -  $^(u>) 


(A- 14) 


where 


U  (w)  =  group  velocity  with  frequency 

^s(u)  =  group  delay  time  at  the  source 

+>?(u)  =  group  delay  time  through  the  instrument 


-A4- 


Since  this  total  group  delay  time  (A-14)  is  frequency  depen¬ 
dent,  energy  at  different  frequencies  arrives  at  different  times 
on  the  seismogram.  Given  that  the  total  energy  of  s(t)  is  equal 
to  that  of  S£w(t)  by  the  same  arguments  presented  earlier,  this 
means  that  the  maximum  energy  density  of  s(t)  will  always  be  less 
than  that  for  sfw(t)  (the  case  where  the  energy  for  every  frequency 
arrives  at  the  same  time).  Thus,  the  signal  is  compressed  in  time 
without  loss  of  energy. 

We  now  show  how  the  noise  is  affected  by  the  matched  filter 
operation.  What  we  compare  is  the  noise 

p  oo 

N(t)  =  j  |N(w)|  exp  +  wt)l  dw  (A-15) 

-w  ^  * 

with  the  matched  filtered  n^ 


ICO 

N(w)  Y*(u>)e1(1)t  dw 

-  -  oo 

or  for  the  whitened  reference  signal  case 

0  CO  „ 

nfw(t)  =  J  I N ( u> )  |  exp[i(-<f0  (w)  +  ♦  4<(w)  +  wt)J  dw 

The  group  delay  at  frequency  w  in  (A-17)  is  given  by 


(A- 16) 


td(")  ‘  +  <>oy(w)  (A* 

Since  1(w)  is  random,  its  derivatives  f(w)  are  also  random. 

The  second  two  terms  on  the  right  side  of  (A- 18)  introduce  a  syste¬ 
matic  shift  in  group  delay  time  with  frequency,  but  because  Y(w)  is 
random,  t^(w)  is  still  random  and  consequently  the  noise  is  not  com¬ 
pressed.  Since  the  matched  filter  operating  on  n(t)  in  the  whitened, 
reference  signal  case  conserves  energy  by  the  earlier  arguments) 
and  is  a  linear  filter  with  a  flat  amplitude  spectrum,  the  noise 


-AS- 


has  the  same  mean  (zero)  and  the  same  variance  (RMS  amplitude)  as 
before,  so  tha;  the  mean  energy  density  in  n£w(t)  is  the  same  as 
the  mean  energy  density  in  n(t).  See  Papoulis,  196S,  p.  34S-347 
for  a  proof  of  this  point  (stationarity  is  assumed). 

Thus  we  conclude  that  the  matched  filter  operation  increases 
the  signal  energy  density  while  not  increasing  the  noise  energy 
density.  Tnis  means  an  enhancement  in  any  conventional  time  domain 
S/N  estimate.  As  discussed  in  the  text  we  can  take  further  advantage 
of  this  effect  in  terms  of  higher  S(w)/N(w)  in  calculations  of 
spectra.  The  gain  results  from  having  essentially  all  the  signal 
power  but  less  total  noise  power  at  all  frequencies  contained 
in  a  time  window  short  compared  to  the  dispersed  signal  duration. 

The  procedure  is  to  compute  the  spectrum  of  a  window  of  length  T 
centered  on  the  peak  of  the  matched  filtered  output.  Mathematically 
this  amounts  to  multiplying  in  the  time  domain  by  a  unit  amplitude 
box-car  function  of  length  T;  the  corresponding  operation  in  the 
frequency  domain  is  a  convolution  with  a  function  of  the  form  s 
This  results  in  a  smoothed  spectrum;  the  distortion  produced  due 
to  smoothing  depends  both  on  T  and  the  shape  of  the  signal  spectrum. 
Figure  4C  illustrates  this  effect  for  the  signal  spectrum  used  in 
the  test  cases  of  this  study.  Generally  speaking  the  signal  band¬ 
width  controls  the  acceptable  minimum  values  of  T  because  the  approxi 
mate  form  of  the  matched  filtered  signal  is  as  can  be  seen 

from  setting  S(uj)  =  1  over  a  finite  frequency  band  2'Af  in  equation 
(A- 1 2 ) .  Other  well  known  types  of  time  domain  windows  could  be 
employed  which  would  produce  different  smoothing  operators,  but 
the  general  effects  on  the  spectra  would  be  similar  to  the  case 
discussed  here. 

Since  the  case  where  the  reference  signal  spectrum  of  the  matchet 
filter  is  not  whitened  is  equivalent  to  prefiltering  x(t)  with  a  phast 
less  filter  whose  spectrum  is  f Y  ( j  ]  and  then  using  the  whitened 
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reference  signal,  all  the  arguments  above  with  respect  to  increasing 
signal  energy  density  relative  to  the  noise  energy  density  hold  for 
this  case  as  well. 

As  a  final  note  it  should  be  pointed  out  that  all  the  analysis 
done  in  this  Appendix  assumes  infinite  limits  in  both  time  and 
frequency,  while  in  practice  we  must  work  with  finite  times  and 
finite  frequency  bands.  However,  it  can  easily  be  shown  that  the 
results  of  this  section  apply  equally  well  to  the  actual  situations 
involving  finite  intervals,  provided  the  proper  care  is  taken  to 
avoid  aliasing  and  wrap-around  in  computing  the  Fourier  transforms. 


APPENDIX  B 


STATISTICAL  ANALYSIS  OF  SPECTRAL  METHODS 


Consider  the  spectral  representation  of  the  operations  involved 
in  the  array  estimates  of  spectra;  namely  if  x^(t)  =  s^(t)  +  n^(t)is 
of  length  Tj,  and  n2i(t)  is  a  noise  sample  of  length  T2  ahead  of 
the  signal  on  the  i*h  channel  then,  denoting  e.g.,  |A(w)|2  by  A2(w), 

T 

[s2(co)]i>n  =  l  2  { [sCco)  +  N^O))]2  -  yi.  N^(u)}i  ( B*  1 ) 

i"l  2 

sum  of  individual 

station  spectra  (Methods  I  and  II) 

^  K  K  T 

[  S2  (w)  ]  1 1 1  =  [y  .2  (S(U))  ♦  N^cu)^2  -  [£  2  (yi  N2(w)).]2  (B-2) 

*  spectrum  from  beamed  trace  (Method  III) 

We  will  show  that  both  approaches  give  valid  estimates  of  the 
signal  spectrum  but  that  the  latter  (Method  III)  is  preferable  because 
it  always  has  a  smaller  variance.  The  presentation  below  follows 
the  analysis  of  this  problem  by  R.  Shumway  (personal  communication). 

To  simplify  the  statistical  analysis  let  us  assume  that  Tj  =  T-, 
on  each  channel;  that  the  true  noise  power  spectrum  is  the  same  for 
all  channels;  that  the  noise  is  uncorrelated  across  channels  and  on 
each  channel  it  is  random  and  stationary;  and  finally  that  the  signal 
spectrum  is  identical  on  all  channels.  Any  of  these  except  the 
assumption  that  the  noise  is  random,  stationary  and  uncorrelated  can 
be  relaxed  without  altering  the  basic  results  concerning  the  perfor¬ 
mance  of  the  two  approaches.  Under  the  conditions  just  stated: 

x  i  ( t )  =  s  ( t )  ♦  nH(t)  i  =  1 ,  K  (B-3) 
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so 


X  i  C^)  *  S  (to)  +  N ^  ( oj ) 


i  =  1 ,  K 


(B-4) 


and 


n2i(t)  -  N2i(oJ) 

where  N.^(u))  and  N2^(w)  are  the  finite  Fourier  transforms  of  each 
of  the  realizations  of  the  random  noise  process. 

E{|  Ni±  (w)  |  2}  =  t||N'2i  (oj)  |  2}  =  N2(o;)  for  all  i  (B-S) 

I  X  i  I  2  =  |  S  (u> )  *  Nli(u))|2  =  S20)  ♦  2  Re  [s*((o)Nli(w)l  (B-6) 

+  .N  £  ±  C  oj  ) 

so 


E{|x(w)  |2}  =  S2  (oj)  +  \'2(uij 


(B-7) 


since 


li{Re  (S*N)}  =  0 

Using  (B7 )  and  (B5)  we  find  that  if  the  noise  is  time  stationary 

E{|Xi|2  -  j  N  ^  i  |  “  }  ^  |  S  (  oj  )  |  2  _  single  station  (B-8) 

signal  power  estimate 
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2 

and  for  an  array  estimate  of  S  (oj)  we  get 


K 

J  2  {lXil2  *  |N2i|2}  «  | S (w) | 2  (Methods  I  and  II  (B-9) 

i  °  1 

The  variance  is  given  by 

K 

var  |  S  (oj  )  |  2  =  -2  f,  [var  |  xi  (<*j)  I  2  ♦  var  |  N2  A  I  2  ]  (B-lO) 

K  i  =  1 

If  we  assume  that  X. (w)  (complex)  is  distributed  normally  then  it 
is  represented  as 


X.  (w)  -  N(S(w),  N  (w) ) 


[Re  S  |  N2/ 2  0  \  1 

Im  S,  \  0  ,\V2  /  J 


(B-ll) 


Then  |X^|  is  distributed  as  chi  square  with  two  degrees  of  freedom 
and  non-centrality  parameter  <52  where  5 2  (to )  *  2 

j  N  ( to )  |  “ 


|xi(u)r  .  x2 

N 2  ( Cl) ) /  2  2’  6 


(B- 12) 


so  that 


|Xi(w) | ‘ 

N  2  (id)  /  ? 


=  4  +  46 4  (u>) 


(B-13) 
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or 


var  |  X±  (go)  !  2  =  N4(w)[l  ♦  S2(gj)]  =  N4(w)[l  ♦  2S2  (u>)  /N2  (w)]  (B-13) 


likewise 


var  | N 2 ^ (w) j  2  =  N4  (“) 


(same  as  above  with  6>»o) 


Thus,  from  these  results  and  (BIO) 

K 


var  | S  (w)  |  j  n  =  I2  2  2N4(u))  [l  ♦  ^1] 
1,11  k  i  =  1  1  IT(ui)  J 


(B-14) 


=  |  N4(w)[l  ♦  S~2~^  1  (Methods  I  and  II) 

in  (to)  J 


Now  consider  Method  III  where 


|S(U>)  I  fu  =  2  x.(w)|2  -  1^2  N2i(«)| 


=  |x(w)r  -  I^TT^JI 


2  (difference  of 

squares  of  means) 


Then 


E{irr£7i2}  •  e{|s(U)  +  2 } 

=  S 2  (to )  +  2E  {r4  +  Efl^T^TI2} 

,  s2m  * 


(B-15) 


(B- 16) 
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and 


e{ittjwi2}  -  Uiifldl 


so 


e{in>iT|2  -  llTjn^TI 2}  =  !S(U))| 


To  find  the  variance  first  note  that 


TOT  ~  N  (S  (to) ,  N2  (to)  /K) 


or  (Re  X,  Im  X)  -N 


fReS  N' 
\lmS,\  0 


/  2K  0 


N2/2K 


so 


ll(id 


N  / 2K 


-  x:  2 


2,v  (to) 


chi  square  with  2  degrees 

of  freedom  and  noncentrality 
2  ,  ^ 

parameter  v  (to) 


where 


v2  (to) 


2K  Sa^ 
N^(to) 
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(B- 17) 


(B- 18) 


(B-19) 


thus 


var 


{^1 


4  +  4  v 


(B- 20) 


or 


var 


|X(u))|2  =  (1  +  v2(w)) 


.  N  ip)  (i  ♦  2k  ) 

K2  N2(w) 


and  similarly 


(B- 21) 


var  |TT7(o>)  |  *  N  4^" 
1  IT 


same  as  above  with 


v  =  o 


so 


var  {|X(w)|2  -  |7T>)|2}  =  var  \J\2  ♦  var  |T?2I 


or 


var 


{|i(u)|2}  -  isl^i  [l  ♦  K  j£lid]  (Method  III) 


(B- 22) 
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Then  we  can  summarize  these  results  in  the  following  table: 

TABLE  B-l 


METHOD  STATISTIC 

I,  II  i  £  {IXjW  2  -  INjjWI2} 


MEAN  VARIANCE 

s2(„)  .  Silail 

K  L  NZ(o))J 


III 


or 


N2i  (u))  |  2  |TCui)  |  2 


lx  X  Na<“)l2 

J  =  1 

liqwi2 


*  KsIm] 

k  N^ui) 


The  variances  can  be  compared  by  computing  the  efficiency 


e 


var  III  1  f 1-K  S2/N21 

var-T7TT  *  *  I  'i7“S  Vn7  J 


s2/n2  ♦  E 
1  V  SZ/N2 
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This  quantity  is  bounded  by  o  <  e  <  1  which  shows  that  it  is  always 
better  to  use  Method  III,  that  is,  beam  before  calculating  spectra. 
We  can  further  note  that  the  performance  depends  on  the  signal  to 
noise  ratio  ^uch  that  for  large  signal  to  noise  ratios  the  methods 
have  almost  equivalent  variances  whereas  for  small  S/N  values 
Method  III  has  a  variance  approaching  ^  smaller  and  becomes  highly 
preferable  to  employ.  As  an  example  take  K  =  4,  10,  20,  and  ®  and 
see  how  the  efficiency  varies  with  S/N.  This  is  shown  in  Figure  18. 
It  is  evident  in  all  cases  that  for  S/N  values  below  about  2 
Method  III  gives  much  better  results  (i.e.  smaller  variance). 

The  result  that  Method  III  always  gives  better  (smaller 
variance)  signal  power  estimates  is  of  great  practical  importance 
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because  only  one  pair  (signal  window,  noise  window)  of  spectral 
computations  is  needed  whereas  K  pairs  of  spectra  are  required  in 
Methods  I  and  II;  using  Method  III  clearly  results  in  a  very  great 
savings  in  computing  time. 
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APPENDIX  C 

COMPUTATIONAL  DETAILS 


Program  COLLAPSE,  whose  basic  structure  and  computational 
procedures  are  shown  in  Figure  19,  was  developed  for  use  in  the 
present  study.  This  program  was  used  to  process  the  actual  events 
discussed  and  the  test  cases  in  which  synthetic  signals  generated 
by  Program  NEWSYN  had  been  imbedded  in  noise,  and  to  calculate  all 
the  spectral  estimates.  Options  are  available  for  matched  filtering 
by  a  large  reference  event  which  can  be  adjusted  for  additional  or 
less  dispersion,  a  synthetic  filter  generated  from  dispersion  data, 
and  a  chirp  filter  transformed  from  velocity  data.  The  synthetic 
signals  themselves  were  used  as  matched  filters  for  the  test  cases. 

Fourier  spectra  were  estimated  for  the  first  3400  seconds  on 
the  test  seismograms,  convolved  with  the  Fourier  spectra  estimated 
from  700  seconds  of  the  synthetic  signal  for  each  station,  and 
inverse  transformed  to  produce  the  matched  filtered  outputs.  The 
RMS  value  of  an  1800  second  noise  sample  on  the  test  seismograms 
was  computed  for  each  station,  and  the  matched  filtered  outputs 
were  weighted  by  1/RMS  and  added  together  to  produce  the  «nm  trace. 
The  spectrum  of  a  300  second  sample  about  the  peak  output  of  the 
sum  trace  was  then  corrected  for  noise  with  the  spectrum  of  an 
1800  second  noise  sample  (Method  III).  At  each  station  power 
spectra  were  estimated  for  a  700  second  signal  sample  on  the  test 
seismograms  and  a  300  second  sample  about  the  peak  of  the  matched 
filtered  outputs  and  were  corrected  for  noise  with  the  power  spectra 
of  the  1800  second  noise  samples  on  the  test  seismograms.  The  program 
was  modified  so  that  the  power  spectra  from  the  test  seismograms 
(Method  I),  the  power  spectra  from  the  matched  filtered  outputs 
(Method  II),  and  the  power  spectra  of  the  synthetic  signals  were 
weighted  by  1/ (noise  RMS)  and  then  summed.  The  1/2  peak  to  peak/RMS 
Signal  to  Noise  Ratios  on  the  matched  filtered  outputs  and  the  sum 
trace  were  computed  and  are  shown  on  Figures  5B,  5C,  6B,  6C,  7B, 

7C,  8B,  and  8C. 

Analogous  procedures  were  used  on  the  actual  events  but  with 
appropriate  signal  windows  on  the  input  traces. 
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PROGRAM  DESCRIPTION 


SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 

DIGITAL  COMPUTING  SECTION 


A.  IDENTIFICATION 

Title:  COLLAPSE 

COOP  Identification:  Z  COLLAPSE 
Category: 

Programmer:  D.B.  Rabenstine  and  J.W.  Lambert 

Date:  July  1969 

B.  PURPOSE 


This  program  processes  time  series  data  from  magnetic  tapes  in 
SDL  subset  format.  The  processes  (each  optional)  which  can  be  performed 
on  the  data  include  (1)  time  shift,  (2)  demagnify,  (3)  band  pass  filter, 
(4)  matched,  phase-equalization  or  chirp  filter,  (5)  weighted  summation, 
and  (6)  signal- to-noise  computation.  There  is  also  an  option  to  do  a 
spectral  analysis  of  signals  both  before  and  after  filtering. 

Matched  filters  may  be  read  from  magnetic  tape  in  subset  format  or 
constructed  from  dispersion  data  from  cards.  Matched  filters  read  from 
tape  may  be  modified  by  data  read  from  cards. 

C.  USAGE 


1.  Operational  Procedure:  This  is  a  Fortran-63  main  program  which 
calls  the  following  subroutines:  SYSLIMIT,  LIBDATE,  PI,  ERASE,  FETCH, 
SPL0T,  FACTER,  C00LER,  SMUTHE,  P0LYE1  C00LBACK,  BPFILT,  PEAKT0PK, 
ABML,  TAPER,  and  PLAT. 

2.  Parameters : 


Card  1  (Repeated  for  each  set,  see  ISUM,  card  10) 


Col. 

Format 

Name 

Description  of  Usage 

1-5 

15 

MTYPE 

=  1,2,3,  or  4,  matched  filter,  read  filter  from 

tape . 


---  11,12,13,  or  14,  read  matched  filer  from  tape 
and  modify  with  dispersion  data  from  cards. 


=  21,22,23,  or  24,  construct  chirp  filter. 


Col.  Format  Name 


P 

f 

£. 


Description  of  Usage 

»  -l,-2,-3,  or  -4  construct  synthetic  filter 
(matched  or  phase  equalize)  from  dispersion 
data  from  cards. 

units  digit  -  1,  "normal"  filter  types  as 
indicated  above,  no  spectra. 

units  digit  -  2,  "normal"  filter  types,  compute 
and  plot  power  spectra  (of  signal  in,  noise  in, 
signal  out;  see  IWHITE  below) . 

units  digit  =  3,  filter  response  (type  as 
indicated  above)  divided  by  noise  power  spectrum, 
no  spectra  output. 

units  digit  =*  4,  filter  response  divided  by 
noise  power  spectrum,  compute  and  plot  power 
spectra,  (see  IWHITE  below) 


6-10 

15 

I WHITE 

*  0  normal  filter  response; 

=  1  white  filter  response 

11-15 

15 

IPOP 

=  1  plot  all  traces;  filters,  inputs,  outputs, 
sum. 

=  0  plot  only  outputs  and  sum. 

=  -1  plot  only  sum 

16-20 

15 

L 

The  number  of  time  points  in  both  the  filter  and 
input  series  is  extended  with  zero's  to  2L. 

21-25 

15 

NA 

Number  across  on  Calcomp  plots  (see  IPOP). 

26-30 

15 

NCODE 

Calcomp  plot  time  scale  increment  in  hundredths 
of  an  inch  per  point. 

31-40 

F10.0 

RANGE 

Calcomp  plot  amplitude  in  hundredths  of  an  inch. 

41-50 

F10 . 0 

DT 

Time  series  sampling  interval  in  seconds  (after 
decimation,  if  any). 

51-55 

15 

IN0RM 

=  -1  summation  (if  any)  unweighted 

=  0  summation  with  weights  equal  to  1.0/noise 
RMS)  . 


=  1  summation  with  weights  equal  to  (signal  peak 
amplitude)/ (noise  RMS)^ 


i 


L  -  3 


Col. 

Format 

Name 

Description  of  Usaqe 

50-55 

15 

IFILT 

=  0  no  band  pass  filtering  , 

>0  number  of  points  to  be  specified  on  bandpass 
filter  response  (must  be  specified  from  D.C.  to 
folding  frequency) 

61-65 

15 

ISM00N 

All  power  spectra  are  smoothed  by  a  running 
average  (boxcar)  of  2xISM00N  +  1  points. 

Card 

2  (Repeated  for  each  set) 

Col. 

Format 

Name 

Description 

1-80 

10A8 

I  DENT 

Arbitrary  identification  for  printer  and  Calcomp 
outputs. 

Card 

(s)  3  (For 

each  set 

if  IFILT  >0) 

Col. 

Format 

Name 

Description 

1-10, 

8F10.0 

AMP (I), FREQ  Points  (4  per  card)  on  response  curve  of 

bandpass  filter.  Must  be  specified  in 
order  of  increasing  frequency  from  D.C. 
to  folding  frequency  (Hz).  Response  value 
at  frequencies  between  those  specified  ar 
linearly  interpolated. 


Card  4  (For  each  case. 

if  MTYPE  =  21-24:  chirp  filter) 

Col,  Format  Name 

Description 

1-10 

F10.0 

DI 

Distance  (epicenter-station)  for  which  filter  is 
to  be  designed. 

11-20 

F10.0 

DO 

ristance  for  closest  station  in  this  set  (=DI  if 
single  case  set) . 

21-30 

F10.0 

CH 

Highest  group  velocity  of  this  filter. 

31-40 

F10.0 

TL 

Longest  period  (i.e.  period  corresponding  to  CH) 
of  this  filter. 

41-50 

F10.0 

CL 

Lowest  velocity  of  this  filter. 

51-60 

F10.0 

TS 

Shortest  period  (i.e.  period  corresponding  to  CL 
of  this  filter. 

61-70 

F10.0 

DT 

Sampling  interval  of  filter  (should  be  same  as 
data  after  decimation). 

Card  5 

(For  each  case, 

if  MTYPE  ■  1-14:  matched  filter  from  subset  tape) 

Col. 

Format 

Name 

Description 

1-5 

15 

ISE1S 

Seismogram  number  of  filter. 

6 

A1 

ISTAR 

*if  operator  is  to  mount  new  tape,  otherwise 
blank. 

7-10 

14 

IC 

Channel  number  (order)  of  filter  seismogram. 

11-15 

15 

IS 

Starting  point  of  filter  relative  to  first  point 
of  on  tape  (before  decimation).  If  IS<0,|IS| 
zero's  are  inserted  in  front  of  filter. 

16-20 

15 

IP 

Number  of  points  in  filter  (after  decimation) . 

21-25 

15 

NDEC 

Decimator;  1<NDEC<5. 

26-30 

F5.1 

SCALEF 

Scale  factor  (counts/u)  applied  to  filter.  (Set 
to  1.0  if  left  blank). 

Card  6 

(For  each  case, 
11-14:  filter 

if  MTYPE  =  -1  to  -4:  synthetic  matched  filter;  or 
from  tape  to  be  modified). 

Col. 

Format 

Name 

Description 

1-5 

15 

KC 

>0,  order  of  polynomial  fit  of  phase  velocity. 

=  0,  use  filter  from  previous  case  (synthetic 

M.F.  only) 

<0,  new  distance,  but  use  polynomial  from 
previous  case. 

6-10 

15 

KI 

>0,  order  of  polynomial  fit  of  instrument  phase 

<^0,  no  instrument  phase  correction 

11-20 

F10.0 

DELTA 

distance  1:  for  filter  modification  =  path  length 
increment;  for  synthetic  filter  =  total  length 
of  single  segment  path  or  first  segment  length 
of  two  part  path. 

21-30 

F10.0 

Cf6 

Highest  velocity  of  interest  in  filter. 

31-40 

110 

KC2 

>0,  order  of  polynomial  fit  of  phase  velocity  for 
second  segment  of  two  part  path. 

41-50 

F10.0 

DELTA2 

If  KC2>0,  =  length  of  second  segment  of  two 

part  path. 


c .  r 


of  two 


Card  (s) 

_7  (If 

KC>0) 

Col. 

Format 

Name 

Description 

1-10, 

11-20 

t  »  * 

8F10.0 

CC(I),I«1, 

KC 

Polynomial  coefficients  of  first  segment 
velocity  as  a  function  of  frequency  (Hz^) . 

phas 

Card  (s) 

8  (If 

KC2>0) 

Col. 

Format 

Name 

Description 

1-10, 

11-20 

•  «  « 

8F10.0 

CC2 (I ) , 1=1, 
KC2 

Polynomial  coefficients  of  second  segment 
velocity  as  a  function  of  frequency  (Hz) 

pha 

Card  (s)  9  (If  KI>0) 

Col. 

Format 

Name 

Description 

1-10, 

11-20 

•  •  • 

8F10.0 

CI(I) ,1=1, IT 

Polynomial  coefficients  of  instrument  phase 
shift  as  a  function  of  frequency  (Hz) . 

Card  10  (For  each  case) 

Col. 

Format 

Name 

Description 

1-5 

15 

ISEIS 

Seismogram  number  of  time  series  to  be 
filtered. 

6 

Al 

ISTAR 

*  if  operator  is  to  mount  new  tape,  otherwise 
blank. 

7-10 

14 

ic 

Channel  number  (order)  of  time  series  to 
filtered . 

be 

11-15 

15 

IS 

Starting  point  of  time  series  relative  to 
first  point  on  tape  (before  decimation) . 
IS<0,  ] IS |  zeros  are  inserted  in  front  of 
time  series. 

If 

16-20 

15 

IP 

Number  of  points  to  be  filtered  (after 
decimation) 

21-25 

15 

NDEC 

Decimator;  1  <_  NDEC  <_  5 

c.c> 


Col. 

Format 

Name 

Description 

26-30 

15 

ISUMT 

-  0,  do  not  put  sum  trace  on  save  tape. 

>0/  put  sum  trace  on  save  tape. 

<0,  read  save  tape  label  from  cards  (11-13), 
write  label  on  save  tape  put  sum  trace  on 
save  tape. 

31-35 

15 

ISAVT 

*  0,  do  not  put  this  filter  output  on  save 
tape . 

<0,  put  this  filter  output  on  save  tape 
<0,  read  save  tape  label  from  cards  (11-13), 
write  label  on  save  tape,  put  this  filter 
output  on  save  tape. 

36-40 

15 

ISUM 

=  0,  do  not  add  this  filter  output  into  sum. 
<0,  clear  sum,  then  add  this  filter  output. 
>0,  add  this  filter  output  to  sum. 

>1,  add  this  filter  output  to  sum, 
this  is  the  last  case  in  set. 

41-45/ 

215 

KN1,KN2 

Noise  window  taken  from  point  KNl  to  point 

KN2  (relative  to  point  IS)  for  computing 
signal-to-noise  ratio  and  also  noise  power 
spectrum. 


51-55, 

215 

KS1, KS2 

Signal  window  in  filter  output  taken  from 
point  KSl  to  point  KS2  (relative  to  point  IS) 
for  computing  signal-to-noise  ratio  and 
signal  power  spectrum. 

61-65, 

215 

KE1, KE2 

Signal  window  in  input  time  series  taken  from 
point  KEl  to  point  KE2  for  computing  input 
signal  power  spectrum. 

71-75 

15 

IHILB 

=  0,  normal  output. 

>0,  Hilbert  transform. 

76-80 

F5.0 

SCALEX 

Scale  factor  (counts/p)  applied  to  time  series 

Card  11  (If  ISUMT  or  ISAVT 

<0,  Cards  11-13  are  save  tape  label) 

Col. 

Format 

Name 

Description 

1-10 

110 

K1 

Seismogram  number 

11-20 

110 

K2 

Number  of  channels 

21-30 

110 

K3 

Number  of  points  per  channel 

31-40 

F10.2 

S4 

Sampling  rate 

c -7 


Col. 

Format 

Name 

Description 

41-50 

F10.2 

T5 

Hour  of  first  point 

51-60 

F10.2 

T6 

Minute  of 

first  point 

61-70 

F10.2 

T7 

Second  of 

first  point 

Card 

12  (If, 

ISUMT  or 

ISAVT  <0) 

Col. 

Format 

Name 

Description 

1-3, 

*  r 

25A3 

ICH(K)  , 

K«1,K2  Channel 

identifiers 

4-6, 


Card 

_13  (If, 

ISUMT  or  ISAVT<  0) 

Col. 

Format 

Name 

Description 

1-72 

9A8 

ID(K),K*1,9  Event 

Identification 

3.  Space  Required: 

4.  Temporary  Space:  Tape  only 

5.  Alarms:  If  problems  are  encountered  retrieving  a  seismogram 
tape,  an  alarm  message  will  be  printed: 

*****  nE  a  XY  **»■-* 

where 

X  »  logical  unit  (1  or  2). 

Y  *  1  if'  seismogram  not  on  tape. 

*  2  channel  number  requested  <1. 

=  3  channel  number  requested  >  number  of  channels  available. 

=*  4  decimator  >5. 

6.  Error  Returns:  See  Alarms 

7.  Error  Stops:  See  Alarms 

8 .  Tape  Mountings : 
a) .  Input: 

Unit  1:  If  MTYPE  =  1-14,  mount  subset  tape  containing 
matched  filter,  else  BY  unit  1. 


Unit  2  s 


Mount  subset  tape  containing  time  series 
to  be  filtered. 


b) .  Output: 

Unit  3:  If  ISUMT  t  0,  or  ISAVT  r*  0,  mount  blank 
tape  for  save  tape,  else  BY  unit  3. 

Unit  4:  If  plots  are  desired,  mount  plot  tape,  else 
BY  unit  4. 

Unit  5;  If  units  digit  of  MTYPE  is  2  or  4,  mount, 
tape  for  printer  plots,  else  By  unit  5. 

d) .  Equivalences: 

/E/7=50/6=51 

9.  Formats : 

a) .  Input: 

Card:  See  Parameters 

Tape:  Units  1  and  2  are  subset  format. 

b)  .  Output: 

Card :  None 

Tape:  Unit  3  is  subset  format  if  user  provides  label 
cards  (see  Parameters). 

Unit  4  is  Calcomp  plot  tape  from  SPLOT. 

Tape:  Unit  5  is  special  format  for  use  interval  to 
program. 

10.  Selective  jumps:  None 

11.  Timing:  About  one  minute/case  for  matched  filtering;  two 

minutes/case  for  spectra. 


